from pylab import *
import tables
import cPickle
from scipy.optimize.optimize import fmin_cg, fmin_bfgs, fmin
from scipy.ndimage import filters

data = cPickle.load(open("mnist-train-images-deskewed.dump"))
data.shape = (60000,225)
data[:,0] = 1.0
classes = cPickle.load(open("mnist-train-labels.dump"))
targets = zeros((len(data),10))
for i in range(len(targets)): targets[i,classes[i]] = 1

def sigmoid(x):
    return 1/(1+exp(-clip(x,-20,20)))
    
def sigmoidp(y):
    return y*(1-y)



class MLP:
    def init(self,n,m,d):
        self.A = randn(n,m)
        self.a = randn(m)
        self.B = randn(m,d)
        self.b = randn(d)
    def forward(self,x):
        y = sigmoid(dot(self.A,x)+self.a)
        z = sigmoid(dot(self.B,y)+self.b)
        return z
    def backward(x,target,eta):
        y = sigmoid(dot(self.A,x)+self.a)
        z = sigmoid(dot(self.B,y)+self.b)
        delta_z = 2*(z-target)*z*(1-z)
        B -= eta * outer(delta_z,y)
        b -= eta * delta_z
        delta_y = dot(delta_z,self.B)*y*(1-y)
        A -= eta * outer(delta_y,x)
        a -= eta * delta_y

class Err:
    def __init__(self,n=10000):
        self.n = n
        self.total = 0.0
        self.count = 0
    def add(self,x):
        l = 1.0/self.n
        self.total = (self.total*(1.0-l)+x*l)
        self.count += 1
    def value(self):
        return self.total

def logreg_gd(data,target,eta=0.1,A=None,a=None,niter=None):
    n,m = data.shape
    n,d = target.shape
    if A is None: A = randn(d,m)
    if a is None: a = randn(d)
    e = Err(); d = Err()
    if niter is None: niter = maximum(3*len(data),1000000)
    for i in range(1,niter):
        j = i%n
        pred = sigmoid(dot(A,data[j].T)+a)
        delta = pred-target[j]
        e.add(argmax(pred)!=argmax(target[j]))
        d.add(sum(abs(pred-target[j])))
        delta = 2*delta*pred*(1-pred)
        l = eta * i**-0.5
        A -= l * outer(delta,data[j])
        a -= l * delta
        if i%100000==0: print i,e.value(),d.value()
    return (e.value(),A,a)



def logpred(data,A):
    return sigmoid(dot(A,data.T)).T
def logloss(data,target,A,verbose=0):
    if A.ndim==1: A = A.reshape(target.shape[1],data.shape[1])
    pred = logpred(data,A)
    loss = sum((logpred(data,A)-target)**2)
    if verbose: print "loss",loss
    return loss
def dlogloss(data,target,A):
    if A.ndim==1: A = A.reshape(target.shape[1],data.shape[1])
    pred = sigmoid(dot(A,data.T)).T
    delta = 2*(pred-target)*pred*(1-pred)
    result = dot(delta.T,data)
    return result.ravel()
def logopt(data,targets,start=None,maxiter=100000):
    n,d = data.shape
    n,c = targets.shape
    data = c_[data,ones(len(data))]
    A = start
    if A is None: A = 0.01*randn(c,d+1)
    def f(x): return logloss(data,targets,x,verbose=1)
    def fprime(x): return dlogloss(data,targets,x)
    result = fmin_bfgs(f,A.ravel(),fprime=fprime,maxiter=maxiter)
    result.shape = (c,d+1)
    return result



x = 0.1*randn(1000)
y = (abs(x)<0.5)
aug = c_[x,x**2]

logopt(aug,array([y]))



def logit(p):
    return log(p/(1.0-p))

def wlstsq(a,b,w):
    assert b.shape==w.shape
    n,m = a.shape
    n_,d = b.shape; assert n_==n
    # weighted least squares: (AT C A) x = AT C b
    # the rows of w are diagonal entries of a norm
    rhs = w*b
    lhs 

def logreg_rw(data,targets,eta=10.0,rtol=1e-3,maxiter=1000):
    ntargets = targets
    last = inf
    for i in range(maxiter):
        A,residues,rank,sv = linalg.lstsq(data,ntargets)
        pred = dot(data,A)
        if i==0: print "lerror",sum((pred-targets)**2)
        spred = sigmoid(pred)
        deltas = (targets-spred)
        ntargets = pred + deltas*eta
        error = sum((spred-targets)**2)
        improvement = (last-error)/error
        print i,error,improvement,"/",mean(A),sqrt(var(A)),"/",mean(deltas),sqrt(var(deltas))
        if i>10 and improvement<rtol: break
        last = error
    return A



def ridge_naive(data,targets,l=0.0):
    n,m = data.shape
    return dot(dot(linalg.inv(dot(data.T,data)+diag(l*ones(m))),data.T),targets)



fdata = array([filters.gaussian_filter(g.reshape(15,15),2.0).ravel() for g in data])


